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Importance sampling of trajectories has proved a uniquely successful strategy for exploring rare 
dynamical behaviors of complex systems in an unbiased way. Carrying out this sampling, however, 
requires an ability to propose changes to dynamical pathways that are substantial, yet sufficiently 
modest to obtain reasonable acceptance rates. Satisfying this requirement becomes very challenging 
in the case of long trajectories, due to the characteristic divergences of chaotic dynamics. Here we 
examine schemes for addressing this problem, which engineer correlation between a trial trajectory 
and its reference path, for instance using artificial forces. Our analysis is facilitated by a modern 
perspective on Markov Chain Monte Carlo sampling, inspired by non-equilibrium statistical mechan¬ 
ics, which clarifies the types of sampling strategies that can scale to long trajectories. Viewed in 
this light, the most promising such strategy guides a trial trajectory by manipulating the sequence 
of random numbers that advance its stochastic time evolution, as done in a handful of existing 
methods. In cases where this “noise guidance” synchronizes trajectories effectively, as the Glauber 
dynamics of a two-dimensional Ising model, we show that efficient path sampling can be achieved 
for even very long trajectories. 


I. INTRODUCTION 

Recent advances in non-equilibrium statistical mechan¬ 
ics have given fresh perspectives on computational pro¬ 
cedures applied to fluctuating molecular systems. The 
Jarzynski work relation and the Grooks fluctuation the¬ 
orem, for instance, provide routes to compute equilib¬ 
rium quantities from non-equilibrium measurements [H- 
1^. Here, we demonstrate that traditional Metropolis- 
Hastings Markov Ghain Monte Garlo (MGMG) can be 
similarly viewed as a procedure to extract equilibrium 
sampling from generically non-equilibrium processes. 
Monte Garlo trial moves drive a system away from the 
steady state distribution, and an entropy production can 
be assigned to these driven transformations. This inter¬ 
pretation provides an elegant way to understand chal¬ 
lenges encountered in MGMG sampling, one that is espe¬ 
cially revealing for MGMG sampling of trajectories. Path 
sampling methods suffer routinely from profound ineffi¬ 
ciency when trajectories of interest become long. From a 
non-equilibrium perspective on MGMG, we provide sim¬ 
ple and quantitative ways to understand the inefficiency. 

Importance sampling of trajectories has enabled stud¬ 
ies of a mwiad of dynamical processes in physics and 
chemistry [7Hl2j|. In particular, reaction rates and mech¬ 
anisms can be found by transition path sampling (TPS), 
which examines the subensemble of trajectories that com¬ 
plete a reaction [I^. The practicality of TPS depends 
intimately on the design of the Monte Garlo (MG) move 
set. Namely, the moves must generate correlated trajec¬ 
tories so that a trial trajectory is likely to exhibit sim¬ 
ilar dynamical behavior as the previously sampled tra¬ 
jectory. Chaotic divergence and microscopic reversibil¬ 
ity of equilibrium dynamics informs the construction of 
two such moves, the so-called “shooting” and “shifting” 
moves [l^- These methods generate correlated trajec¬ 
tories by propagating alternative histories from highly 
correlated initial configurations. For sufficiently short 


trajectories, the imposed correlation at one time serves 
to strongly correlate the trajectories at all times. Long 
trajectories, however, are problematic: trial trajectories 
either lose all useful correlation with the reference path, 
or else they coincide so closely with the reference that 
changes are impractically small [l^ . In both cases the ef¬ 
ficiencies of shooting and shifting moves plummet as tra¬ 
jectories grow longer. Sampling trajectories that involve 
slow molecular rearrangements and diffusive processes 
stands to benefit significantly from alternative methods 
of generating correlated trajectories. 

We consider three different ways to guide long tra¬ 
jectories: introducing auxiliary forces; selecting among 
series of short trial segments, as in Steered Transition 
Path Sampling (STePS) [l^; and advancing stochastic 
integrators with correlated random numbers (which we 
refer to as “noise guidance”) [il[i3- Of the three, only 
noise guidance yields an MG entropy production which is 
subextensive in the trajectory length. The other schemes, 
which accumulate extensive entropy production, cannot 
efficiently extend to sampling of long trajectories. Strong 
noise guidance is not, however, a panacea; correlated 
noises need not imply correlated trajectories. We il¬ 
lustrate this point by considering Glauber dynamics of 
a two-dimensional Ising model and Langevin dynamics 
of a two-dimensional Weeks-Ghandler-Andersen (WCA) 
fluid. Only when microscopic degrees of freedom have 
a small number of discrete possibilities, as in the lat¬ 
tice dynamics, is it possible to generate correlated long- 
timescale trajectories by tuning the noise. 

The structure of the paper is as follows. First we 
introduce and discuss the perspective of MG moves as 
non-equilibrium processes which produce entropy, detail¬ 
ing the consequences of constraints analogous to fluc¬ 
tuation theorems and the second law of thermodynam¬ 
ics. Next we review transition path sampling in stochas¬ 
tic dynamics and demonstrate the challenge posed by 
long trajectories in the context of trajectory sampling 
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of a one-dimensional random walker. We then ana¬ 
lyze alternative strategies to correlate long trajectories 
of a one-dimensional single particle system and of a 
two-dimensional Ising model. Finally, we explore how 
strongly noise guidance correlates trajectories in exam¬ 
ple systems, and then conclude. 


II. MARKOV CHAIN MONTE CARLO 
ENTROPY PRODUCTION 


We start by discussing a very general perspective on 
traditional Metropolis-Hastings MCMC sampling 0 
0- Consider the problem of sampling a configuration, x, 
according to probability distribution f’(x). For example, 
X could be a vector of the coordinates and momenta of 
N hard spheres, the state of spins in an Ising model, or 
all coordinates of a classical trajectory. The Metropolis- 
Hastings algorithm generates a Markov chain, which can 
be thought of as a dynamics through configuration space 
with the steady-state distribution P(x). This dynamics 
obeys detailed balance but is not necessarily physical. 

One typically splits each Monte Carlo move into two 
steps. First, a change from x to a new state x is pro¬ 
posed according to a generation probability, Pgen [x —>■ x]. 
Throughout this paper we will refer to x as a reference 
and X as the trial. This trial is conditionally accepted 
with probability 


Paccept[x x] = min [l, e , 


where 


,, , 1„ f (x)P,.„[x ^ it] 

P(x)Pgen[x -)■ x] 


( 1 ) 

( 2 ) 


P{uj). Following the more general demonstration of an 
entropy production fluctuation theorem [2l|, note that 

P(a;) = J dx dx P(x)Pgen(x —?> x)(5(a; — a;(x, x)) 

= J dx dx e“^’‘'’^^P(x)Pgen(x —?► x)S{uj + uj{x, x)) 

= e“P(-w), (3) 

with a;(x, x) representing the entropy produced by a pro¬ 
posal move from x to x and S denoting the Dirac delta 
function. We are more likely to propose moves with 
positive entropy production than we are to choose their 
negative counterparts. The straightforward corollary, 
(w) > 0, is by analogy a statement of the second law, 
and the equality is satisfied if and only if P(u;) = (5(a;). 
A further consequence of Eq. [3] relates the MC accep¬ 
tance rate to the probability of attempting a move with 
a negative value of w, which we call p<. Specifically, 

(Paccept) = y [l’ = 2p<, (4) 

which has been noted in the related context of replica 
exchange Monte Carlo [? ]. As (ui) increases, p<, and 
therefore (Paccept), tends to decrease. We will see that 
(w) scales with the number of driven degrees of freedom, 
such that Monte Carlo sampling of chain molecules or of 
long trajectories becomes especially challenging. 

We focus below on the sampling of dynamical pathways 
(rather than individual configurations). In this case lj is 
an “entropy production” only by analogy, since the “non¬ 
equilibrium transformations” effected by MC trial moves 
occur in the more abstract space of trajectories. Lessons 
from Eq. ([3]) are nevertheless illuminating in the context 
of this abstract space. 


Together, these two steps ensure detailed balance, guar¬ 
anteeing that the equilibrium distribution P(x) is sta¬ 
tionary under the MC protocol. Lacking the conditional 
acceptance step, such an MC procedure would generally 
drive a system away from its equilibrium distribution. 
We find it instructive to view this notional process as a 
genuine non-equilibrium transformation, one that would 
generate nonzero entropy in most cases. In the formalism 
of stochastic thermodynamics, the resulting entropy pro¬ 
duction corresponds precisely to the quantity uj defined 
in Eq. (ED [13. 

The MC acceptance step effectively filters realizations 
of this non-equilibrium process, with a bias towards low 
values of w. By construction, the bias exactly negates 
the tendency of trial move generation to drive a system 
out of equilibrium. From this perspective, the counter¬ 
vailing tendencies of proposal and acceptance are akin to 
the operation of a Maxwellian demon, which by contrast 
filters realizations of equilibrium dynamics with a bias 
that creates a non-equilibrium state. 

Viewing the procedure in the language of entropy pro¬ 
duction distributions reveals an important asymmetry of 


III. TRANSITION PATH SAMPLING WITH 
STOCHASTIC DYNAMICS 

A. Trajectory Space and Trajectory Subensembles 

Let us now specialize to the sampling of discrete¬ 
time stochastic trajectories with a fixed number of steps, 
fobs- The probability of observing a trajectory, x{t) = 
{xq , xi,..., xj^^^}, can be written as 

^obs 1 

Po[x(0] OC Pinit(xo) (5) 

where pinit is a distribution for the initial time point, 
frequently an equilibrium or steady state distribution. 
The probability of each time propagation step is denoted 
p{xt —>■ xj+i), the form of which depends on details of 
the stochastic dynamics. We refer to this propagation as 
the natural dynamics. Representative trajectories can be 
generated by sampling the initial state and propagating 
natural dynamics. 
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In many contexts, it is useful to study a biased 
trajectory ensemble constructed to highlight particular 
rare events. Common examples include the reactive 
subensemble, 

C’reactive[x(t)] OC Pq [x(t)]/lA (xq j/lB ^ (6) 

and the so-called tilted ensemble, 

Ptiited[x(t),s] oc Po[x(t)]e-^^W‘)l. (7) 

In the former case, and he are indicator functions 
which constrain the endpoints of the trajectory to fall 
in regions of phase space corresponding to reactants and 
products of a chemical reaction or other barrier cross¬ 
ing process [l^. In the latter case, iC[x(<)] is an order 
parameter reporting on dynamical properties of the tra¬ 
jectory (e.g. the current [^, activity [l2;[23, or entropy 
production [23 - f^ l and s sets the strength of bias 
These biased ensembles highlight classes of trajectories 
only rarely sampled by the natural dynamics. To effec¬ 
tively sample them, a Markov chain of correlated tra¬ 
jectories is constructed. The correlations between sub¬ 
sequent steps of the Markov chain ensure that newly- 
generated trajectories are likely to share the rare features 
that made the prior trajectory a good representative of 
the biased ensemble. 


B. Sampling with Shooting Moves 

One of the most general and effective methods for gen¬ 
erating a trial trajectory is the shooting move, which is 
particularly well-suited to sampling equilibrium dynam¬ 
ics [l^. The move proceeds as follows. A discrete time, 
4hoot, between 0 and fobs is uniformly selected and des¬ 
ignated the shooting time. The state of the system at 
4hoot, perhaps slightly modified, is then propagated for¬ 
ward and backward in time with natural dynamics to 
yield a trial trajectory, x(f). The probability of generat¬ 
ing this trial takes the form 

.fgen[x(t)] C>Cpinit(xo)_Pgen ^ ^tshoot) ^ 

^shoot 1 ^obs 1 

n pi^t+i -t xt) p{xt xj+i), 

t—O t —^shoot 

( 8 ) 

where p is the transition probability for time-reversed 
dynamics and Pgen is the probability of the perturbation 
at the shooting time. In the language of Section [HI the 
entropy produced by this trial move is given by 

^ jOinit (xq ) /j-A (Xq ) fee )pgen (Xf shoot ^ ^^shoot ) 

Pinit (xo ) hA (xo ) hB (Xi.bs )Pgen (x* shoot ^ ^^shoot ) 

f-' p(xt+l ^ Xt)p(Xt Xj+l) 


For long trial trajectories to be accepted by the MCMC 
scheme, w must be small. However, when p and p are 
not equal (as is the case in driven processes), the sum 
in Eq. has order tobs nonvanishing terms [28l| . Con¬ 
sequently, (w) scales linearly with fobs, and P(w) adopts 
the long-time form 

P{uj) ~ exp [-tobs7(w/tobs)], (10) 

with large deviation rate function /(w/tobs)- From this 
asymptotic expression for the entropy produced by a TPS 
move, one might generally expect that the correspond¬ 
ing acceptance rate decreases exponentially as fobs grows 
long. 

This extensive growth of (w) with time has an impor¬ 
tant and general exception, namely the case of microscop¬ 
ically reversible dynamics. Under those conditions, the 
sum in Eq. ([5]) vanishes and the only entropy production 
is contributed from the endpoints of the trajectory (e.g., 
/lA and /ib). Since this entropy production is subexten- 
sive in time, long trajectories appear no more difficult 
to sample than short ones. The acceptability of trial 
trajectories, however, is also subject to biases like those 
expressed in Eqs. m and ©• Because long trajectories 
typically decorrelate strongly from one another, the rare, 
biased qualities of a reference trajectory (e.g., reactivity 
or inactivity) are recapitulated in the trial path with a 
probability that also decays with tobs- 

We conclude that the challenges for efficiently sampling 
long trajectories are twofold. The TPS move must pro¬ 
duce entropy that is subextensive in observation time or 
the method will not scale to long trajectories. Addition¬ 
ally, one must preserve strong correlations between x(t) 
and x(t), so that rare properties of interest are retained 
in the trial trajectory. In the next section we show that 
these two goals are often conflicting. In particular, we ex¬ 
amine three general schemes for engineering correlations 
between reference and trial trajectories in shooting-like 
moves. Two of the schemes fail to exhibit subextensive 
entropy production scaling while the remaining scheme 
can only maintain strong trajectory correlations in spe¬ 
cial cases. 


IV. GUIDED DYNAMICS OF A ID RANDOM 
WALKER 

We explore the three methods for trial trajectory 
generation in the specific context of a one-dimensional 
discrete-time random walker with equation of motion 

Xt+l=Xt+^t, (11) 

where at the timestep, the noise is drawn from 
the normal distribution with zero mean and variance cr^. 
As a simple illustration focusing on the effects of entropy 
production, suppose we want to sample the unbiased tra- 
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FIG. 1. Three guidance schemes for generating a trial trajectory that maintains proximity to a reference trajectory. For 
the specific case of a one-dimensional random walker, upper panels illustrate the consequences of (a) artificial corralling forces 
(b) preferential selection of short trial branches and (c) correlated noise histories. Bottom panels show the corresponding 
distributions of trajectory space entropy production u). The intensity of red shading reflects the probability that trial moves 
are rejected. For cases (a) and (b), the average entropy production is nonzero and grows with trajectory length tohs- With 
an appropriately designed noise guidance scheme (c), symmetric selection of noise variables results in identically zero entropy 
production for all trajectory lengths. 


jectory distribution 


A. Guiding Forces 


Po[x{t)] oc (5(a;o)exp 


E 




{xt+l - XtY 

2a2 


( 12 ) 


where the initial position is set to zero without loss of 
generality [1^. To construct a reference trajectory a;(<), 
we draw a value for S,t at each timestep and propagate the 
walker’s position according to Eq. (HU. A trial trajectory 
is then generated by evolving dynamics from the same 
initial configuration (with a different realization of the 
noise or perhaps even a different equation of motion). 

We imagine that it is desirable for the trial trajec¬ 
tory to retain a significant correlation with the reference 
path. This goal is motivated by the challenges of sam¬ 
pling biased ensembles as discussed above, but for the 
sake of simplicity we do not include such a bias here. To 
ensure this correlation, we employ shooting moves that 
differ from the conventional procedure described in Sec¬ 
tion UTIB] Specifically, we implement and scrutinize three 
distinct ways to engineer correlation over long times: (a) 
adding artificial forces that pull the trial trajectory closer 
to the reference, (b) preferentially selecting among sets 
of otherwise unbiased short path segments, or (c) using 
correlated histories of noises. We assess the influence of 
these three biasing methods on the MCMC efficiency by 
characterizing the distribution P{uj). 


We first consider effecting correlations with guiding 
forces, i.e., artificial contributions to the effective poten¬ 
tial that tend to lead the trial trajectory toward the refer¬ 
ence. This strategy is equivalent to using steered molecu¬ 
lar dynamics to generate new trajectories. The trial 
trajectory x(t) is grown with the equation of motion 

Xt+I = Xt + it + k{xt - Xt). (13) 


We denote as the trial trajectory noise at timestep t, 
also drawn from a Gaussian with mean zero and variance 
a^. The linear spring constant k adjusts the strength of 
correlation between reference and trial trajectories. The 
probability that this guided dynamics generates a partic¬ 
ular trial from the reference is given by 


Pgen[x{t) — >■ i(t)] OC exp 


E 




(it +1 - Xt - k{xt - Xt)f 
2a2 


(14) 

The entropy production associated with the trial move 
depends also on the probability of generating the reverse 
TPS move, growing the reference trajectory with extra 
forces pulling it close to the trial. It is straightforward 
to compute uj from Eq. 


^obs 1 

W =-y {xt - Xt){Xt+l + Xt+l - Xt - Xt). (15) 
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In this approach, w can be physically interpreted as the 
difference between two work values: that expended by 
the artificial force to guide the trial trajectory, versus 
the work that would be required to conversely guide the 
reference. The resulting distribution of entropy produc¬ 
tion, obtained from numerical sampling, is shown in in 
Fig. (Ha). 

Since uj is given by a sum over all tobs timesteps, 
P{u}) adopts a large deviation form as in Eq. (ITUl) . and 
(w) (X tobs- These scaling properties are demonstrated 
numerically in Fig.[2][a) and analytically in Appendix [A] 
In the appendix we re-express uj in terms of the ^ and ^ 
variables, which can be integrated over to yield 

[(2 - k) -l + {k- . (16) 


Indeed, for 0 < k < 2, this expression gives the antici¬ 
pated long time scaling with tobs, 


(w) 


2A:tobs 

2-k' 


(17) 


As seen in Fig. [TJa), the negative-w tail of which 

gives rise to MCMC acceptances, becomes correspond¬ 
ingly small for large tobs- 


B. Guiding Choices 

In both Sections IIII Bl and IIV Al we showed that time- 
extensive entropy production arises generically when we 
do not use natural (forward) dynamics to generate a tra¬ 
jectory. Dynamical biases can alternatively be achieved 
by preferentially selecting among different examples of 
natural dynamics. At a high level, conventional TPS [l^ 
is just such an approach, constructing biased trajectory 
ensembles through selection rather than artificial forces. 
Can this strategy be used effectively to impose resem¬ 
blance between reference and trial trajectories? 

We consider a scheme very similar in spirit to the 
STePS algorithm fla. L ike configurational-bias MC sam¬ 
pling of a polymer [^, the STePS procedure generates 
a long trajectory by piecing together short segments, as 
illustrated in Fig. [TJb). To generate a new segment, one 
starts at the end of the previous segment and samples 
a collection of short, unbiased trajectories according to 
the natural dynamics, which we will refer to as branches. 
One of these branches is selected as the next segment of 
the trial trajectory, with a preference for branches that 
stay close to the reference trajectory. (Proximity could 
be judged in different ways, e.g., through Euclidean dis¬ 
tance in the full phase space, or with respect to an order 
parameter). Though each branch is grown with natural 
dynamics, the added segment is biased. To show that 
this bias affects acceptance rates in the same manner as 
the guiding forces bias, we compute the entropy produced 
by a TPS move. 


Starting at the initial condition of the reference trajec¬ 
tory, we grow n branches of length r according to 

(18) 

where a is an index over the n independent samples of 
the natural dynamics. Of these n possibilities for the 
segment of the trial trajectory, we select branch a with 
probability 


Tselect(tr) — 


fi\4r 


(“)| 






(19) 


where / is a weighting function with a maximum when 
its argument is zero, for example a Gaussian centered on 
zero. The reference trajectory is indicated by a super¬ 
script (0). Starting from the end of the chosen branch, 
the growth procedure is repeated with n new branches of 
length r. 

While each time propagation step uses segments of 
unbiased natural dynamics, the selection of preferred 
branches exerts a bias which ultimately leads to a non¬ 
vanishing entropy production. 


^obs/”^ 

w = - ^ In 
2 = 1 




( 20 ) 


where ai is the index of the selected branch for the 
segment. The calculation of this entropy production re¬ 
quires generation of the backwards TPS move, in which 
the (0) branch is always selected. 

In the preceding section we discussed that the entropy 
produced by guiding forces could be thought of in terms 
of a work performed by the bias. From that perspective, 
this guiding choices scheme trades work for information. 
We bias the dynamics not by applying explicit forces but 
instead by selecting the preferred branches based on in¬ 
formation about the likelihood that a branch stays close 
to the reference. In particular, a; is a difference between 
the Shannon information associated with selecting the set 
of trial branches which produced the trajectory x(t) and 
the information associated with selecting the reference 
branches in a reverse TPS move. As with biasing forces, 
the trajectory space entropy production exhibits large 
deviation scaling with (w) oc fobs- Numerical demon¬ 
strations of this scaling are provided in Fig. HKb). Con¬ 
sequently, acceptance probabilities drop precipitously in 
the long time limit. 


C. Guiding Noises 

As a third scheme for engineering path correlation, we 
consider generating a trial trajectory with natural dy¬ 
namics but with biased noises. Rather than trying to 
corral trajectories to proceed along similar paths, one 
may impose much simpler correlations between their un¬ 
derlying noise histories [l^ . Consider a TPS trial move 
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^/^obs ^/^obs 


FIG. 2. Entropy production statistics for trial moves with guiding forces (a) and guiding choices (b), as discussed in 
Sections IIV Al and IIV Bl of the main text. Results in (a) are shown for fc = 0.1; black curves indicate the long-time behavior 
determined analytically in Appendix 1X1 Results in (b) are shown for n = 3, r = 10 and with f{x) = In both panels, 

different colors indicate different trajectory lengths, ascending from left to right: fobs = 30 (red), 50 (orange), 100 (yellow), 200 
(green), 500 (cyan), and 1000 (purple). 


which consists of re-propagating dynamics from the ini¬ 
tial timestep using new noises ^ that differ only slightly 
from the old noises 

6 = + \/l - (21) 

where rji is sampled from a Gaussian distribution with 
zero mean and variance a^. In Section (IIVBI) . the sym¬ 
bol a was used as an index, but here we redefined a to 
be the parameter controlling noise correlations. Guiding 
Gaussian noise variables in this manner has been referred 
to as a Brownian tube proposal move 0. Unlike the 
prior two kinds of moves, the Brownian tube proposal 
produces strictly vanishing entropy production ui for all 
trials regardless of trajectory length. The cancellation 
results from some algebra after writing the path weights 
and generation probabilities in terms of the noise vari¬ 
ables. 




exp 

1 

2^ 

2(l-a2)o.2 

P(6Pgen(|^6 exp 

e ■ 

exp 

(€-«€)= ■ 

20-2 

2(l-a2)CT2 


( 22 ) 

where ^ is a vector detailing noises at all times. 

The vanishing entropy production is achieved by in¬ 
dependently sampling the noise variables. In the previ¬ 
ously discussed approaches, the bias applied to any one 
noise variable depended on how far astray the trial tra¬ 
jectory had drifted from its reference up to that point in 
time. Such history-dependent biasing coupled the sam¬ 
pling of one noise variable to all of the previous noises, 
ultimately giving rise to the nonvanishing uj. By sam¬ 
pling all noises independently, we can perturb the ^ vari¬ 


ables in a symmetric manner. For noises drawn from a 
Gaussian distribution, this symmetric perturbation was 
given in Eq. (ED, but the strategy of symmetrically sam¬ 
pling independent noises generalizes to other choices of 
stochastic dynamics. For example, Hartmann has ap¬ 
plied these methods with uniform random variable noises 
to Monte Garlo dynamics in the form of Wolff dynamics 
of a two-dimensional Ising model [32| . 

Using correlated noise histories to sample Monte Garlo 
trajectories avoids the time-extensive bias that arose 
from guiding paths in configuration space, but this merit 
comes at a cost. If the reactive or tilted ensembles are 
to be sampled, it is important that the guidance scheme 
produces highly correlated trajectories. That is to say, 
the X coordinates, not just the ^ coordinates, must be cor¬ 
related. When will similar noise histories produce similar 
trajectories? In the remainder of the paper we address 
this question in the context of two dynamical systems, 
one on-lattice and the other off-lattice. 


V. EFFICACY OF NOISE GUIDANCE 

In the preceding section we noted that sampling tra¬ 
jectories with noise-guided shooting moves avoids a time- 
extensive MG entropy production. However, we seek cor¬ 
related trajectories, not just correlated noises. When tra¬ 
jectories with correlated noises synchronize, efficient path 
sampling of long trajectories can be achieved. But un¬ 
der what conditions should such synchronization be ex¬ 
pected? We investigate this question by studying lattice 
dynamics of a two-dimensional Ising model and off-lattice 
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FIG. 3 . Correlations between a reference trajectory and a trial trajectory generated by the noise guidance method described 
in Section fV Al using eair = Csite = 10 “® and Cacc = 0 . 1 . The two trajectories begin with identical initial conditions and evolve 
through 100 sweeps of push up/push down Monte Carlo steps at pJ = 0 . 3 , h = 0. Final configurations of reference and trial 
trajectories are shown in (a) and (b), respectively. The site-wise overlap between these two configurations is depicted in (c), 
where black indicates spin alignment and white indicates anti-alignment. 


dynamics of a WCA fluid, also in two dimensions. We 
show that synchronization can be achieved with a suit¬ 
able treatment of Ising dynamics. This success does not 
extend to our example of off-lattice dynamics. 

A. Ising Dynamics 

Let us first consider a two-dimensional Ising model 
consisting of N spins. The I**' spin, denoted ai, takes 
the value ±1. The lattice evolves, at inverse temperature 
P, under single spin-flip Glauber dynamics with Hamil¬ 
tonian 

H = -h (23) 

* (b) 

The spins interact in the usual Ising manner; they cou¬ 
ple to nearest neighbors with coupling constant J and 
to an external field h. Each spin-flip trial move requires 
us to choose two random numbers uniformly from [0,1). 
One random number, ^site, determines which site will 
be flipped. The other random number, ^acc, determines 
whether to accept or reject the flip. Given ^site and ^acc, 
the spin-flip move is deterministic: 

1. Ghoose spin i = ceiling(^siteA^) to act on. 

2. Gonstruct a trial state by flipping spin i. 

3. Compute the energy difference, AE, between the 
original configuration and the trial. 

4. Accept the spin flip if ^acc < (1 + exp(/?AiiJ))“^. 

By carrying out fobs sequential MG sweeps, each consist¬ 
ing of N spin flip moves, we construct an Ising trajectory, 
a{t). The effective unit of time is thus taken to be a MC 
sweep. 


Now consider a noise-guided trial TPS move designed 
to generate a trajectory d{t) which is correlated with 
<j{t). At every MC step we alter ^site and .^acc to some 
trial values, ^site and ^acc- There is significant freedom in 
doing so, while producing zero entropy [3^ . and we ana¬ 
lyze one particular choice. We focus first on updating the 
noise that chooses which spin to flip. With probability 
1 — Csite we reuse the old noise, i.e., ^site = Csite- Other¬ 
wise, we uniformly draw a new value of ^site from the unit 
interval. The tunable parameter Csite controls the corre¬ 
lation between noise histories of the reference and trial 
trajectory. We update the noises that control conditional 
acceptance, ^acc, in an analogous manner. Another pa¬ 
rameter, Eacc, is the probability of drawing new noise for 

Cacc- 

Starting with the initial configuration of cr(t), we con¬ 
struct a{t) by performing spin flips with the new trial 
noise history. The trial and reference trajectories start 
in identical configurations, but we expect the correlation 
to decay as MG time advances. To monitor the similarity 
between reference and trial, we study the site-wise prod¬ 
uct between cr and a as illustrated in Fig.|3l The average 
of this product over all spins, 

1 ^ 

CT • CT = — ^ o^di (24) 

' i=l 

is a measure of correlation between a and d. Decorre- 
lated configurations return a value of zero while identical 
configurations return one. Fig.|4l(a) shows that the corre¬ 
lation between a{t) and a(t) decays to zero at long times. 
The rate of this decay is tuned by Cacc and Csite, the pa¬ 
rameters controlling the extent of noise correlation. At 
long times (cr • a) eventually approaches zero, even with 
the strongest noise guidance. The corresponding “uncor¬ 
related” configurations, however, bear a subtler resem- 
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(a) (b) 




FIG. 4. Average overlap between reference and trial trajectories of a 40 x 40 two-dimensional Ising model with /?J = 0.3. 
Results are shown in (a) for ordinary Glauber spin flip dynamics, and in (b) for the modified directional dynamics described in 
Section FV Al Two trajectories with identical initial conditions and different-but-correlated noise histories (solid lines) maintain 
a nonzero steady state overlap at long times only for the case of push up/push down dynamics. The same steady state values 
are obtained when the two trajectories evolve from very different initial conditions (dashed lines) generated by independently 
assigning each spin at random. Different colors indicate different values of the noise guidance parameter tacc. Ensemble-averaged 
results are shown for Csite = Cdir = 0.001, with averages performed over 500 independent pairs of trajectories. 


blance. Some regions are significantly correlated, while 
others are significantly anti-correlated, averaging to give 
a ■ a m 0. 

Motivated by this subtler resemblance, we introduce a 
minor alteration in the implementation of Glauber spin- 
flip dynamics. Below we detail this modification and 
show that it does in fact enable the preservation of cor¬ 
relation between trajectories over very long times. As 
observed in the context of damage spreading, different 
choices of Ising model Monte Carlo dynamics can result 
in identical equilibrium states yet different dynamical 
properties . The effectiveness of correlated noises 

in guiding trajectories is, in effect, one such dynamical 
property. 

In particular, we replace step 2 of the spin-flip move to 
include a directionality. We introduce another random 
number, ^dir G [0,1), used to decide the trial state. If 
^dir < 0.5, the trial state is a down spin; otherwise it is up. 
Whereas the conventional trial move effects an attempted 
spin flip, this trial move can be viewed as an attempt to 
push the spin either up or down, depending on the state 
of ^dir- As with the other noise variables, correlations 
between ^dir and ^dir are tuned by the probability Cdir of 
resampling the noise. 

The addition of directionality to the spin-flip dynam¬ 
ics results in moves which are trivially ineffectual. For 
example, half of the spin moves attempt to “push up” 
a spin which is already up. These moot moves are ab¬ 
sent from the traditional implementation of single-spin- 
fiip Glauber dynamics, which attempts a spin flip at ev¬ 
ery step of MC time. In every other respect, the two 


schemes generate Markov chains with identical statistics. 
They can therefore be made identical by excising moot 
moves, or, equivalently on average, by scaling time by a 
factor of two. 

As Fig. HKb) illustrates, the push up/push down im¬ 
plementation of single-spin-flip Glauber dynamics allows 
the trial trajectories, a{t) to remain tunably close to a{t) 
for long times. By incorporating information about spin 
change directionality into the noise history, the noises 
signal not just how likely a spin is to change, but in what 
direction it will change. Appropriately chosen e parame¬ 
ters can create trajectories which remain tunably close to 
each other for arbitrarily long times. When averaged over 
the whole lattice, steady state correlation is maintained, 
but the correlations are not spatially homogeneous. As 
MG time progresses, the regions in which two trajecto¬ 
ries are highly correlated move throughout the lattice, 
ensuring ergodic exploration of the trajectory space. 

For push up/push down dynamics, noise guidance does 
not merely preserve correlations that existed at time zero. 
We find that correlated noise histories can in fact induce 
synchronization between trajectories. To illustrate this 
synchronization effect, we have characterized the corre¬ 
lation between trajectories that share similar noise histo¬ 
ries only intermittently. As shown in Fig. [SJ such paths 
acquire similarity during periods of strong noise guid¬ 
ance. This similarity degrades during periods without 
noise guidance, but can be recovered by re-introducing 
guidance, regardless of how significantly correlations have 
decayed. Indeed, even very different initial configura¬ 
tions, propagated with correlated noises, become more 
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— eacc=10“^ 

— eacc=10“^ 

— eacc=10“^ 

— eacc=10° 


FIG. 5. Correlation between reference and trial trajectories of a 40 x 40 Ising model with /3J = 0.3. Plotted lines are averages 
over 500 independent pairs of trajectories that evolve by push up/push down dynamics. Noise histories of each pair are generated 
in a correlated manner with Cgite = Cdir = 10~® during the intervals t = 25 — 50 and t = 150 — 200 (the shaded regions); noise 
guidance is absent at all other times. The site-wise correlation between an example trajectory pair (with Cacc = 10“®) is shown 
above the graph, with time advancing from left to right and adjacent configurations separated by 25 MC sweeps. 


similar with time, their ensemble-averaged correlation 
{a ■ a) converging to the same value as for trajectories 
that are identical at time zero. 

A nonzero steady state value of (cr • a) is the quanti¬ 
tative signature of synchronization. The origins of this 
finite asymptotic correlation are transparent in the limit 
of weak coupling, /3J = 0. With the additional simplifi¬ 
cation h = 0, each attempted spin flip is accepted with 
probability 1/2 based on the value of ^acc, regardless of 
the states of neighboring spins. In this case the steady 
state overlap can be calculated analytically. To do so, we 
derive an equation of motion for the probability p{t) that 
a given spin has identical values in the reference and trial 
trajectories after r MC steps. Note that r differs from 
time t by a factor of N. The long-time, steady-state 
limit of this time evolution, pss = limT-_>ooP(T) , yields 
(ct • o')gg = 2pss — 1. In Appendix [b 1 we tabulate the 
various ways that a selected spin can become identical 
in reference and trial trajectories after a single timestep. 
From this enumeration, and the corresponding probabil¬ 
ities, we find 

(25) 

The various terms of Eq. [51] describe the different ways 
that a single MC move can impact the state of an ar¬ 


bitrarily chosen spin in trial and reference trajectories. 
Since each move of our MC dynamics acts on a single 
site of the lattice, some moves do not involve the tagged 
spin at all, but instead some other lattice site; the hrst 
term in Eq. [55] reflects this possibility. The second term 
accounts for the decrease in overlap when reference and 
trial trajectories accept a spin-flip at the same site but in 
opposite directions. The third term results from the con¬ 
structive action of correlated noises on the tagged spin, 
either maintaining existing correlation or inducing syn¬ 
chronization, as detailed in Appendix |B] The final term 
accounts for random alignment of the tagged spin de¬ 
spite uncorrelated noise variables, a possibility particular 
to degrees of freedom with a limited number of discrete 
states. Equating p(r) and p(r -|- 1) gives the steady state 
probability. 


Pss — 


1 _ (1 _ (1 _ ) 
2 - (1 - ^) (1 - Esite + 


(26) 


Analytically calculating steady-state overlap at finite 
temperature is not straightforward. Numerical results, 
shown for f]J = 0.4 in Fig. |6l indicate that the depen¬ 
dence of overlap on strengths of noise perturbation is 
generically similar to the PJ = 0 case analyzed above. In¬ 
creasing j3J from zero does, however, slow the rate of con¬ 
vergence to the steady state, while decreasing the degree 
of steady state correlation. For all coupling strengths we 
have examined, (cr • can be made arbitrarily close to 
unity by decreasing the various e parameters. This level 
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FIG. 6. Steady state correlations between reference and trial trajectories of a 40 x 40 Ising model, as a function of noise 
guidance strength. Results are shown for push up/push down dynamics with eair = 0.001. The high-temperature limit Eg. 1261 
is shown in (a). Finite temperature behavior (b) was obtained by sampling 4 x lO'^ pairs of trajectories for 340 MC sweeps 
each, with /3J = 0.4. 


of control ensures that one can generate trial trajecto¬ 
ries which are correlated with a reference for all times, 
an essential capability for efficient path sampling of long 
trajectories. 


B. WCA Dynamics 


The success we have achieved in synchronizing Ising 
dynamics with noise guidance should not be expected for 
complex dynamical systems in general. We demonstrate 
this limitation for the specific case of a two-dimensional 
WCA fluid evolving by underdamped Langevin dy¬ 
namics. The purely repulsive particles are propagated 
using an integration scheme that requir es g enerating a 
collection of Gaussian random variables . These 

noises are guided by a Brownian tube proposal, Eq. 

The similarity between trial and reference noise histories 
is controlled by a parameter a that ranges from zero (no 
noise guidance) to one (complete noise guidance). 

Starting from identical initial configurations, we prop¬ 
agate dynamics with correlated noise histories and mon¬ 
itor the difference between trial and reference as 


(|x-i|) 



2=1 


(27) 


where and are the positions of particle i in the 
reference and trial and j-j is the two-dimensional Eu¬ 
clidean distance. At short times, the difference between 


trial and reference trajectories is small, but this differ¬ 
ence grows exponentially, a hallmark of chaotic dynamics. 
Even with exceptionally strong noise guidance, trajecto¬ 
ries cannot be held arbitrarily close to each other for long 
times, as shown in Fig. [71 

Why are we unable to guide the evolution of WCA 
particles as effectively as we guided Ising dynamics? A 
principal difference between the two systems is the like¬ 
lihood of spontaneous local recurrences. In either case, 
trajectories with similar initial conditions but different 
noise histories wander away from one another in a global 
sense, eventually exploring very different regions of con¬ 
figuration space. Correlating their noise histories gen¬ 
erally acts to defer, but not defeat, this divergence. In 
the case of Ising dynamics, however, a given small block 
of adjacent spins will occasionally align spontaneously in 
two trajectories with reasonable probability. Following 
such a spontaneous, local recurrence, correlated noises 
again work to hold trajectories close. Global spontaneous 
recurrence of a large number of spins is, of course, highly 
improbable, but the noise guidance seems to “lock in” 
local correlations every time they spontaneously reoccur. 


VI. CONCLUSION 

Transition path sampling has proven useful for a vari¬ 
ety of equilibrium, as well as non-equilibrium, problems 
in chemical dynamics. The problem of sampling long tra¬ 
jectories, particularly those with multiple intermediates, 
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FIG. 7. Divergence between reference and trial trajectories 
of a two-dimensional WCA fluid with Brownian tube noise 
guidance of strength a. Average distance between the two 
trajectories, as defined in Eq. 1271 is shown for a system of 400 
particles with mass m and diameter a, in a square box with 
side length 24 cr (i.e., density pa^ — 0.694). Underdamped 
Langevin dynamics was propagated with inverse temperature 
/3 = 0.2 and friction coefficient 7 = 0.1 using a timestep 
of 0.005 T, where r = and e is the Lennard-Jones 

interaction energy scale [ 3 ^ . Data are averaged over 500 in¬ 
dependent trial trajectories. 


has hindered a variety of extensions and applications of 
the methodology. We have outlined a modern physical 
perspective from which to assess and address these chal¬ 
lenges. We have demonstrated and discussed successful 
trajectory guidance in the case of Monte Carlo dynam¬ 
ics of an Ising model. Substantial difficulties remain for 
systems with continuous degrees of freedom. 

Our results suggest that effective noise guidance of long 
trajectories requires a nonnegligible probability of spon¬ 
taneous local recurrence, i.e., a significant chance that 
reference and trial trajectories transiently align within 
small regions of space. Such synchronization could be 
particularly helpful for sampling reactive trajectories 
that traverse metastable intermediate states, for example 
the coarsening or assembly of colloidal systems as they 
organize on progressively larger scales. In such cases, 
trial trajectories in the course of path sampling should 
maintain correlations with the reference while passing 
through the intermediates, not just at the endpoints. 
Even without identifying metastable configurations, cor¬ 
related noises could be applied during some intervals but 
not others. A tendency to synchronize would enable trial 
trajectories to explore widely during unguided periods, 
but to be reined in globally by intermittent guidance. 

We anticipate that these noise-guidance methods will 
be effective for other lattice systems as well, but their 
usefulness could depend sensitively on the exact manner 
in which the noise influences dynamics. In particular, 
without incorporating directionality into proposed spin 


changes, we were not able to guide long Ising trajectories. 
Furthermore, Ising dynamics can exhibit spontaneous re¬ 
currence, i.e., transient local alignment between two tra¬ 
jectories regardless of their noise histories. Because small 
blocks of Ising spins can adopt only a modest number of 
configurations, such random local synchronization occurs 
with an appreciable probability. The probability of recur¬ 
rence will likely be lower for models with a larger collec¬ 
tion of possible local configurations, e.g., a Potts model 
or an Ising model with more neighbors. We thus expect 
that the application of noise-guided path sampling could 
face substantial challenges for long trajectories of these 
more intricate lattice models. 
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Appendix A: Entropy Production Statistics for a 
One-Dimensional Random Walker with Guiding 
Forces 

1. Mean Entropy Production 

Here we analytically characterize the entropy produc¬ 
tion distribution, P{uj), for shooting moves generated 
with guiding forces as discussed in Section IIV Al It is 
useful to first rewrite Eq. © in terms of the noise vari¬ 
ables, ^ and For a one-dimensional random walk, the 
position Xt+i and the difference between reference and 
trial trajectory, Xt+i — Xt+i^ can be compactly expressed 
in terms of the noises. 

t 

Xt+l = ^ ^ (-^ 1 ) 

u—0 
t 

Xt+i-Xt+i='^{l-kY~'^{^u-iu)- (A2) 

u—0 

After straightforward algebra it is possible to express uj 
as 

- = (A3) 
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where for convenience we have defined = ^t+^t + 2/i, 
C = it-it, and 


t-i 

(A4) 

u—0 


The main text presents results for a random walk with¬ 
out drift, i.e., with ^ drawn from a distribution with 
mean zero. Here we consider the more general case 
with nonzero mean /r. Noting that {ii = {it~iu) = 
2a‘^5tu and {it iu) = Oj the average entropy production 
is found to be 


tobs-i t-i 

(w)=2 ^ (A5) 

t—0 u—0 I 


The two geometric series are summed to yield 

(w) = 1 + (^ “ i)2tobs] ^ (A6) 

When /c > 2, (w) grows exponentially in fobs- This su- 
perlinear scaling results from coupling between trajecto¬ 
ries so strong that the trial trajectory rapidly tends to 
infinity due to a numerical instability, much like the in¬ 
stability that arises in conventional molecular dynamics 
simulations performed with an excessively large integra¬ 
tion timestep. For the useful range of coupling strength, 
k < 2, (uj) oc fobs in the long time limit. The marginal 
k = 2 case is well-behaved ((w) = 4tobs (^obs — l))j but 
uninteresting for our purposes. 

2. Cumulant Generating Function 

The behavior of the higher-order cumulants can be ex¬ 
tracted from the cumulant generating function In 
This average requires integration over all of the Gaussian 
^ and ^ variables at all times, which can be performed 
inductively. We define (/)(A, /, g, h) as 




f 1 


2t 


exp 


\ 2ay/iT 

■f (-{iU?-{iYi? 


dit ■■■ dit-idi^ ... exp 




2=0 


1 ^^(-{itf-{iT? 


+ xsYt - ^sf 


+ A {St.Yt -1 - hSl, + 2gk{l - 


(A7) 


The integral (j) is defined such that (e“'^“) = 

^(A, 1, 0,1, fobs)- By introducing f,g, and h we can de¬ 
rive recursion relations as we sequentially integrate out 
Gaussian noises at the latest remaining timestep. In par¬ 
ticular, integration over i^_i then it-i returns an inte¬ 
gral of the same form. That is to say 4>[X, fi,gi,hi,t) = 
4>{\, fi+i,gi+i,hi+i,t - 1) with 


fi+i — 


h 


gi+i — X — hi ~\- 


4Agf(l — k^k'^ 


1 — AXgik'^ 


(A8) 

(A9) 


h,+i = i-{i-kY {{x-h,) + 


4Ag^(l — k^k" 


1 — AXgik^ 


(AlO) 


where go = 0 and 


gi+i — A — 1 -l- 


(1 - kYgi 

1 — AXgik'^ 


(A12) 


The scaled cumulant generating function is then given by 
lim — In/e-^"^) =-iln(l-4A5*fe2), (A13) 

^obs ^OO ^Q]3g 2 

where g* is a fixed point of the map given in Eq. (IA12I) . 
Specifically it is the lesser of the two roots of the 
quadratic equation obtained when gi = = g* is in¬ 

serted into Eq. (IA12I) . The numerical Legendre transform 
of this scaled cumulant generating function gives the solid 
black curve in Fig. [5t^a), toward which the results of nu¬ 
merical sampling should converge for long fobs- 


Iterating the map fobs times corresponds to integrating 
over all of the 2tobs integrals in Eq. dm. After some 
algebraic simplification, 

ln(e-^-) = -l 5] ln(l - 4A5 ./c 2), (All) 

i=0 


Appendix B: Ising Model Steady State Correlations 

To derive the steady state correlation between refer¬ 
ence and trial trajectories, we examine the time evolution 
of the probability p{t) that reference and trial overlap at 
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site i and MC step r. Without loss of generality, we focus 
on a particular site, i = 1. Push up/push down moves 
can be grouped into four classes: (i) the reference and 
trial each act on spin 1, (ii) the reference acts on spin 1 
while the trial acts on a different spin, (iii) the trial acts 
on spin 1 while the reference acts on a different spin, or 
(iv) neither reference nor trial acts on spin 1. For each 
case, we derive a transition matrix which maps the vec¬ 
tor (p(r), 1 — p{t)) to its state at MC step r -I- 1. The 
full transition matrix for a step of dynamics is the sum 
of these transition matrices, weighted by the probability 
of each case, 

(Bl) 

The transition matrix for case (iv) is the identity matrix, 
since this case cannot alter the overlap at site 1. The 
transition matrices for cases (i), (ii), and (iii) are Q,R-, 
and S, respectively, the forms of which we now derive. 

When reference and trial act on different spins, only 
one copy of spin 1 (the reference or the trial) can change 
its state. When the two copies differ at site i after r 
steps, overlap is induced with probability 1/4 (i.e., the 
probability that any given move results in a change of 
spin state). For initially aligned copies, loss of overlap 

I 


similarly occurs with probability 1/4. This logic applies 
equally well to cases (ii) and (iii), so 

0 = « = (v( 3/() ■ <“> 

When both reference and trial act on site 1, we must 
account for correlated influence on the two copies. As a 
result, S depends on Cacc and Cdir- To enumerate these 
correlated changes, we denote states of spin 1 at step 
T in the reference and trial as cti and di, respectively. 
After the MC step, these spins are given by and d). 
Table U lists the possible transformations which result in 
overlapping spins (cr( = d') after t + 1 steps. Collecting 
terms in the table and making use of the fact that S' is a 
probability-conserving transition matrix, we find 

(B3) 

Propagation according to the transition matrix T gives 
the overlap probability after a single MC step: 


f p{t+ 1) \ ( Pi^) \ 

yi-p(T-ti)y \i-p{T)j 


(B4) 


The first row of this matrix equation reads, after some 
algebra. 


p{t -f 1) = p(t) - -f 


^site 


N 


1 ^site 


N 


2N 


2N 


1 - 


ir /-I ^acc \ , /-■ ^acc \ / \ 

r(‘-—j+d-—jpM 


(B5) 


We are interested in the steady state solution, found by by N, followed by algebraic simplification, yields Eq. [55] 
settingp(T) =p(r-|-l). Multiplying the equation through of the main text. 
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-1111 
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1-10 1 

-1111 

(^) 
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-110-1 


1 -11-1 

-1 -1 * -1 
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